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ABSTRACT 

In this paper we revisit the problem of the tidal interaction occuring between a proto- 
stellar accretion disc and a secondary point mass following a parabolic trajectory. 
We model the disc response analytically and we compare our results with three- 
dimensional SPH simulations. 

Inviscid as well as viscous hydrodynamics is considered. We show that in a viscous 
system the response derived from inviscid considerations is predominant even for the 
highest estimates of an anomalous disc shear viscosity. The angular momentum lost 
from the disc during the encounter is derived from linear theory, for distant fly-bys, 
as well as the changes to the disc orientation expected in non-coplanar encounters. 
It is shown that the target discs can become warped and precess by a small amount 
during non-coplanar encounters. This small precession is shown to give rise to a relative 
tilt of the disc which is always more important for determining its final orientation 
than is the change to the orbital inclination. 

We discuss the implications of our results for protostellar accretion discs and planetary 
systems. 
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1 INTRODUCTION 

Low mass stars are known to form in compact groups within 
. 5^ I the dense cores of giant molecular clouds. Recent obser- 
■ vations of nearby star-forming regions have revealed that 
$_( ' young stellar objects (YSOs) have a binary frequency in ex- 
cess of that found for field stars and roughly half of the 
current sample are associated with optically thick, geomet- 
rically thin, circumstellar accretion discs (see recent reviews 
by Mathieu 1994, Lin & Papaloizou 1996 and references 
therein). If the mean free path in a star-forming core is suf- 
ficiently small so as to allow close encounters between YSOs 
on a timescale shorter than the lifetime of a protostellar disc 
then the tidal interaction between a secondary object and 
the disc of the primary may be a dynamically significant 
event in the disc's history. 

The possibility that YSOs might tidally interact through 
their circumstellar discs has lead to recent speculation that 
binaries may form as the products of a capture process in- 
volving star-disc encounters (Larson 1990). However it was 
quickly shown that even for the most compact star-forming 
environments known, disc penetrating encounters leading 
to capture would be rare (Clarke & Pringle 1991), with 
the probability of such an encounter being approximately 
one tenth over a disc lifetime. Also observational data im- 
plies that most YSO binaries have components of similar 
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ages (Hartigan, Strom & Strom 1994, Brandner & Zinnecker 
1997), which supports local formation of the components 
rather than capture. We note however that more distant en- 
counters are near certainties for protostellar accretion discs. 

These kinds of star-disc interactions have been extensively 
studied in the context of star formation. Heller (1993) inves- 
tigated numerically with SPH the rate of disc tilting result- 
ing from non-coplanar disc penetrating encounters. Clarke 
& Pringle (1993) investigated numerically, but with a sticky 
particle code, the affect of a disc penetrating encounter on 
the distribution of disc material. Ostriker (1994) applied 
linear perturbation theory to an inviscid hydrodynamical 
model, using a spherical harmonic expansion of the tidal 
potential for distant encounters, to derive an asymptotic ex- 
pression for the angular momentum lost from the disc in the 
encounter. Korykansky & Papaloizou (1995) used a Fourier 
expansion in azimuthal modes, for the tidal potential, and 
calculated the angular momentum exchange without asymp- 
totic assumptions. Hall, Clarke & Pringle (1996) carried out 
a numerical investigation of the disc response using a re- 
duced three-body method with non-interacting particles. 

Due to the complexity of the tidal problem the analytical 
results presently available in the literature, although exten- 
sive, are also unweildy and the numerical studies are chiefly 
qualitative. In this paper we present analytical and numeri- 
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cal (SPH) calculations of the disc response to a parabolic en- 
counter. We present a simplified analysis isolating the most 
significant parts of the response which have been identified 
in the works listed above. We obtain explicit expressions for 
the angular momentum exhange (and other properties of 
the response) and we test these results directly against the 
simulations. 

In Section 2 we give some of the equations basic to our sub- 
sequent considerations. In Section 3 we review concepts of 
the linear disc response and perform a linear mode analysis 
of the fluid equations. In Section 4 we derive the linear disc 
response for distant coplanar encounters and in Section 5 
we extend this to consideration of the response for distant 
non-coplanar encounters. Our numerical method is outlined 
in Section 6 and in Section 7 we present our numerical re- 
sults. In Section 8 we summarise our findings and discuss 
our results within the context of recent observations. 



the midplane value being given by Cg — Q'^H^/2n. 
Following a thin disc approximation we shall also consider 
the disc surface density E, defined by: 

The constant C„ is given by the beta function /3(i,n + 1). 
We do not consider viscous forces explicitly in the equation 
of motion, but where required we make use of the standard 
Q-prescription due to Shakura & Sunyaev (1973). In this 
prescription the usual coefficient of kinematic shear viscos- 
ity, i/, is written as = aCsH. In this case a is the dimen- 
sionless parameter to be determined. The viscous evolution 
timescale is ~ r^/i^, being much larger than the dynamical 
timescale, which is also the timescale on which vertical 
hydrostatic equilibrium can be established. 



2.2 The secondary 



2 BASIC EQUATIONS 

2.1 The accretion disc model 

We consider a non-self-gravitating gas-dynamical model for 
the accretion disc. Processes that may give rise to internal 
heat generation or transport are not considered. Instead we 
adopt a simple polytropic relationship between gas pressure, 
P, and density, p, giving a constitutive equation of state: 



Kp 



l + l/r. 



(1) 



where n is the polytropic index and K is the polytropic 
constant. The associated barotropic sound speed, Cs, is given 
by 
2 dP 

dp 

The constancy of K in this model requires that we assume 
any dissipated energy to be lost from the system. This is 
equivalent to assuming an efficient cooling mechanism oper- 
ating in the disc. 



2.1.1 Equilibrium structure 

Referring to a set of cylindrical coordinates (r, 0, z) with 
the 2-axis coincident with the rotation axis of the disc, the 
standard equation of vertical hydrostatic equilibrium in the 
thin disc approximation is: 



1 9P 



(2) 



where Q — f2(r) is the Keplerian angular velocity of disc ma- 
terial in an initially axisymmetric disc. For a disc satisfying 
the polytropic equation of state (jl|) integration of equation 
(|) gives 



p{r,z) 



2K{n + l) 



{l-z^/H^)", 



where H — H{r) is the total vertical semi-thickness for 
which p{z = iiH) — 0. We may now compute the sound 
speed in the disc in terms of the midplane value, Cs'. 



Cs = Cs(l - Z^/H'^), 



Initially we shall consider the case in which the disc mid- 
plane and the orbital plane of the secondary coincide. Our 
results are then extended to include the case in which the 
two planes are mis-aligned. We take the primary and the sec- 
ondary to be point masses. Neglecting any energy exchange 
phenomena that might occur between the secondary and the 
disc the energy per unit mass of the secondary referred to 
the non-inertial frame centred on the primary is written: 



E = — 
2 



(J- G(A/p + Ms) . o- 



+ 



(3) 



In a coplanar configuration the position of the secondary at 
time t is D{t) = [cit), 9{t),0\. Mp is the mass of the primary, 
Ms is the mass of the secondary and we denote with a dot the 
derivative with respect to time. For a parabolic trajectory we 
set E = and integrate (^) with respect to time. Referring 
to a set of Cartesian coordinates {x,y, z) we can choose the 
y-axis to coincide with the longitude of pericentre. In this 
case the position of the secondary as a function of time is 
given parametrically: 

£) = g[4p,(l- V),01, 



in which 



sinh 



■ sinh 



The distance of closest approach of the secondary to the 
primary is denoted by q and the magnitude of the angular 
velocity of the secondary at pericentre by ujo: 



_ y/2G{Mp + Ms) 

^2 ■ 

The natural origin of time corresponds to the instant of peri- 
centre passage, i.e. D = |il>| — q at t — 0. 



2.3 Tidal potential for distant encounters 

We consider a geometrically thin disc, initially in a state of 
vertical hydrostatic equilibrium and radial centrifugal equi- 
librium, governed by the central potential vl'o alone. This 
initial state is then perturbed from equilibrium by the tidal 
forcing due to a secondary point mass encountering the pri- 
mary/disc on a parabolic prograde orbit. The part of the 
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total potential that is due to the secondary, 'i'' , is consid- 
ered to drive small perturbations of the disc's equilibrium 
structure. Within the context of a fluid model we can then 
apply a linear perturbation analysis of the fluid equations. 
This will always be valid if the affect of the tide is sufficiently 
weak, and amounts to considering the distance of closest ap- 
proach q to be large when compared to the outer radius of 
the disc _Ro- The total potential is 'I' = \I'o -I- where the 
potential due to the primary acting at a point with position 
vector r is given by 



*o = 



GMp 



and the potential due to the secondary, referred to the origin 
which is fixed to the primary, is given by 

GMs GMsr ■ D 



r D 



(4) 



Expanding up to second order in terms proportional to r/D, 
and ignoring terms proportional to (z/r)^, we may then 
write the part of the total potential due to the secondary 
as (Papaloizou & Terquem f 995) 

/■I .-> I T' - # y 

1 " ^:t^ + 



D 



3(r • D) 



2D2 



(5) 



3 THE DISC RESPONSE 

3.1 Wave generation and propagation 

For tightly wound spiral density waves propagating with a 
radial wave number of magnitude \kr\ in a geometrically thin 
self-gravitating disc the standard dispersion relation (Lin & 
Shu f964) is written: 



27tGE|A:J + \krfct. 



(6) 



In the non-self-gravitating limit this may be expressed as 

\kr\ = -Vs^-l, (7) 



where s (Binney & Tremaine 1987) gives the ratio of the 
forcing frequency in a frame rotating with pattern speed Sip 
to the natural radial frequency k (also called the epicyclic 
frequency) : 



s ■■ 
and 

2 _ 
K = 



m{Q — f^p) 



1 f 4o2\ 



The positive integer m is the azimuthal mode number. Res- 
onances occur at s = 0, —1 corresponding to the corota- 
tion resonance and the inner and outer Lindblad resonances 
respectively. When a fluid disc is subject to a perturbing 
force which is stationary in a frame rotating with a con- 
stant angular frequency Qp, a stationary disturbance is set 
up in that frame, the disturbance having the same order of 
azimuthal symmetry as the perturbing force. 
The radial wave number given by (^) indicates short wave- 
length spiral density waves which may be either inward or 
outward propagating. Furthermore, standard results from 
the analysis of the solutions of the dispersion relation rtd) are 



that for non-self-gravitating discs there exists an evanescent 
zone about the corotation radius which occupies the region 
between the inner and outer Lindblad resonances. Inside the 
evanescent zone wave propagation is formally disallowed. We 
can further deduce that the spiral density waves become in- 
creasingly tightly wound as they travel away from the Lind- 
blad resonances. 

In the tidal problem there is a torque exerted on the disc by 
the secondary. The torque principally interacts with the disc 
through the Lindblad resonances (Goldreich & Tremaine 
1979), provided that they are either in, or lie sufficiently near 
to the disc (Lin & Papaloizou 1993) . The applied torque ex- 
cites short trailing waves at the Lindblad resonances which 
propagate away from the evanescent zone with a wavelength 
oc H. As long as the waves remain linear then they propagate 
with a conserved wave action. Inward propagating waves 
then become increasingly tightly wound and increase in am- 
plitude until non-linear damping takes place and the angular 
momentum they transport is deposited in the disc. This is a 
well known process that may lead to induced accretion near 
the disc centre. The net angular momentum transport due 
to the corotation resonance is much smaller and depends on 
the ratio 4f2fc^/E (Korykansky & Papaloizou 1995), evalu- 
ated at the resonant location. 

Consideration of the full problem for parabolic fly-bys re- 
quires analysis of the disc response to an infinite spectrum 
of forcing frequencies. Recent work in this area (Ostriker 
1994, Korykansky & Papaloizou 1995) has shown that a sin- 
gle component of the response predominates: the part of the 
tidal potential with m — 2 azimuthal symmetry applied at 
pericentre. 

Assuming that the secondary-disc interaction occurs as an 
impulse near pericentre we shall treat the disc response 
as resulting from a non-resonant interaction. In this time- 
independent calculation the secondary is assumed to be fixed 
at pericentre so that the torque exerted on the disc can be 
calculated in a simple way. In practice the interaction has a 
finite width in the sense that a resonant perturbation may 
be set up at some time before pericentre passage. As a result 
a phase lag between the density response of the disc and the 
forcing potential is set up, allowing a non- vanishing net tidal 
torque to be applied to the disc. 

We shall determine the disc response with the secondary 
fixed at its pericentre position using a linear mode analysis 
of the fluid equations and apply the calculated torque over 
the estimated time interval of the interaction. Korykansky & 
Papaloizou (1995) report that the interaction interval they 
observe in their numerical models is about one orbital pe- 
riod at the outer edge of the disc, centred at the instant of 
pericentre passage. 



3.2 Linear mode analysis 

The vertically averaged inviscid fluid equations (see Pa- 
paloizou & Lin 1995b, for a recent review) yield radial and 
azimuthal components of the momentum equation: 



dt ^ dr r dS 



VI 



9t 9r r 9(A r 



j_9P _ 9* 
E 9r 9r ' 

1 9P 19* 

Er h<h r dS ' 
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and the continuity equation: 

^ + 19(E.K) + l^(Erl/.)=0. 
ot r or 0(j) 

Note that within the context of the vertical averaging proce- 
dure physical quantities are formally replaced by their ver- 
tically averaged analogs. The radial and azimuthal velocity 
components are denoted by Vr and respectively. 
We proceed by making a standard linear analysis of the ver- 
tically averaged fluid equations such that quantities are per- 
turbed by a small amount from their equilibrium states (e.g. 
E — > Eo -I- S' in standard notation). Then subtracting the 
equations for equilibrium and neglecting quantities of sec- 
ond and higher orders in the perturbation one obtains the 
fluid equations linear in perturbed quantities: 



9K 
dt 

9^ 

dt 

and 

9S' 
'dt 



2nvl 



h * 

9r I Eo 



or 



1 9 / P' ^ , 
h * 

r 9«i I Eo 



r or am r om 



These equations describe the linear response of the model. 
We analyse for azimuthal modes such that we make replace- 
ments like E' Ea(r) exp(mi(/>), noting that physically 
meaningful results are obtained by taking the real parts. 
The complex amplitude functions for the time-independent 
forced velocity components are (Goldreich & Tremaine 
1979): 



GMe 



1-l-if ^ ) (1 -3cos2(?!)) 



The axisymmetric parts of this yield no net tidal effect (Pa- 
paloizou & Pringle 1977). The non-axisymmetric part indi- 
cates a disturbance with azimuthal mode number m = 2. 
Carrying out the potential expansion to terms of order 
{r/qf' yields terms describing modes with m = 1 and m = 3, 
indicating that the dominant response comes from the m = 2 
component with the m = 1 and m = 3 modes becom- 
ing increasingly significant in closer encounters. This is in 
agreement with the findings of Korykansky & Papaloizou 
(1995), those authors claim a satisfactory fit to their numer- 
ical results can be calculated by considering only the Fourier 
components of the perturbative potential corresponding to 
modes m = 1, 2 and 3; the contribution from m = 1 and 3 
modes being much smaller than that from the m = 2 mode. 
Considering a potential perturbation 

* = - — 7T-r exp(2i(/)), 
4 

the velocity and surface density perturbation amplitudes are 
then given by: 



5 GM, 



and 

Ea = 



4 ng3 



(13 - 4n) mEq 3 
: T:—r ■ 



and 



1 

A 



A 



20, 9r 



2mn ^ 9 

h mil — 

r or 



where A = ~ m^Q,^ . Additionally we assume a thin disc 
(i.e. H/r is a small quantity) for the duration of the en- 
counter so that we can neglect the pressure perturbation in 
comparison to the perturbation to the potential. The per- 
turbation to the surface density is then determined from 
the perturbed continuity equation. In terms of the velocity 
perturbations, 



Here Eo = Eo(Po/r)"-' 
Eo = C„7?o 27f(„ + i) 



in which 

" jj^ 2n+l 



In the above; fi denotes the mass ratio Ms/Mp, the rota- 
tional velocity at the outer edge of the disc is Qo ~ 0{Ro) 
and we have considered H/r to be a constant (in practice 
H/r is a weak power law in r, Larwood & Papaloizou 1997). 
We can now deduce the part of the net torque acting on the 
disc which gives rise to a change in the magnitude of the 
disc angular momentum J. The rate of change of the disc's 
angular momentum is then given by (Papaloizou & Pringle 
1977): 



1 

mil 



.Eol^,;,a 19, 

mi h - — (EorK, 

r r or 



We note that the denominator. A, can vanish for m = 1 
azimuthal modes, since for strictly Keplerian rotation n = Q. 
This leads to a singularity in the disc response equation 
which is usually dealt with by employing a dissipative phase 
shift in the density response (see below). 



4 LINEAR DISC RESPONSE FOR DISTANT 
COPLANAR ENCOUNTERS 

Evaluating ^ at pericentre with the plane of the sec- 
ondary's orbit coplanar with the disc midplane we have 



dJ 

dt 



E'-^-rdrd<j!>. 



Notice that the linear terms in the perturbed torque are 
identically zero and that the net torque is second order in 
perturbed quantities. Angular momentum transfer will then 
be a correspondingly small effect in inviscid discs. In fact this 
part of the torque will itself only be non- vanishing if the per- 
turbation to the density becomes phase shifted with respect 
to the perturbing potential such that E' oc expmi(0 — rj). 
In order to obtain an order of magnitude estimate of the 
angular momentum exchange we consider the torque due to 
the leading non-axisymmetric term of the potential expan- 
sion outlined above. Integrating over the disc for the m = 2 
response, and from t = —nQ^^ to t — 7tf2^^, we find 
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AJ; 



J -Jo 
Jo 



37t(13-4n)(7-2n)^2 



16 



8- 



— 1 sin2r;. 



(8) 



Note that we identify the angular momentum of the unper- 
turbed disc, Jo, with the angular momentum content of an 
axisymmetric polytrope: 



Jo = In 



Eor r^rdr. 



In general a density perturbation is set up in phase with 
the potential perturbation at a finite time before pericentre 
passage. It is the differential rotation between the pattern of 
the perturbation at the outside of the disc and the secondary 
that results in a phase shift between the perturbing poten- 
tial and the density response. The existence of a phase lag 
allows a net angular momentum transfer to occur between 
the secondary and the disc. 

Our viscous numerical models also require consideration of 
the angular momentum exchange due to viscous dissipation. 
This has been derived by Papaloizou & Pringle (1977): 



dJ 
dt 



So£ 



rdrd(^. 



(9) 



where $ = 2v{eijeij — ^e?,) defines the rate of dissipation 
of mechanical energy per unit mass in terms of the rate-of- 
strain tensor, eij, of standard fluid theory. Evaluating 
in the simple case of constant viscosity, and for the m — 2 
perturbations calculated above, gives a contribution [AJ]^ 
to the angular momentum exchange: 



[A,/]. 



-707t. 



7 - 2n 
'15-271 



Comparison of this result with equation implies that for 
reasonable values of the polytropic index viscous torques be- 
come important when r] ^ 107?."^. Where Ti. is the Reynolds' 
number (= r^Q.jv) evaluated at the outer edge of the disc. 
Since viscous timescales longer than ~ lO^yr are inferred in 
protostellar accretion discs (Beckwith et al. 1990) we would 
require a dynamical phase shift that is necessarily small com- 
pared with unity. There is no reason why we should expect 
such a small value for hence in general the viscous con- 
tribution to the process of angular momentum exchange is 
expected to be negligible, as is conflrmed by the numerical 
results to be presented below. 



5 EXTENSION TO DISTANT 

NON-COPLANAR ENCOUNTERS 

Assuming the disc to be axisymmetric and the secondary- 
disc interaction to occur only at pericentre, we need only 
consider the inclination angle 5 of the vector pointing to 
pericentre with respect to the initial disc midplane. So if 
we choose the a;-axis to lie along the line of nodes, then TD 
is given by a simple rotation of the previously considered 
coplanar trajectory about that axis: 

D = g[4p, (1 - 4p^)cos5, (1 - 4p^)sin5]. 



The perturbative potential may then be recalculated and 
decomposed into a linear sum of contributions corresponding 
to odd terms in z (denoted by ■i/'o) and even terms in z 
(denoted by t/ie): 

Now, to terms of first order in 2, equation (^) leads us to 
consider non-axisymmetric components 

3 GMs 2 2 r /o- J.\ 

~ 5— r COS exp(2i(p) 



and 



.3GM, 



L 25 exp(ii 



The non-axisymmetric part of the non-coplanar response 
that is even in z, t/ip, carries an additional factor of cos'^ 5 in 
comparison with the analogous coplanar term. Hence inclin- 
ing the orbit diminishes the strength of the response found 
in the coplanar case, being minimised in orthogonal configu- 
rations. Thus the effectiveness of tidal truncation of the disc 
and the strength of non-axisymmetric waves is reduced. In 
addition our estimate for the angular momentum exchange 
is modified: A J — > AJcos^ &. 

The part of the potential that is odd in z, ^o, yields a dis- 
turbance with azimuthal mode number m = 1; this henA- 
mg mode is found to warp the disc and cause its preces- 
sion in circular orbit calculations (Papaloizou & Terquem 
1995, Larwood et al. 1996). But note that the amplitude we 
compute above is two times larger than that for the analo- 
gous zero-frequency term in a circular orbit calculation with 
D — q (see Papaloizou & Terquem 1995). Note also that we 
have recovered a warping potential proportional to sin 25, 
Ostriker (1994) showed that this result holds asymptotically 
for distant encounters. 



5.1 Precession, warping and changes to the 
inclination 

We consider only the warp due to the part of the perturba- 
tive potential with m — 1, i.e. the longest wavelength mode. 
This warp will be the longest-lived and most significant as 
long as the disc response is linear. Also this component will 
cause precessional effects. 

As before we calculate the forced response assuming the 
perturbing potential to be fixed at its pericentre value. We 
suppose that our coordinate system precesses at a rate Up, 
which corresponds to the precession of the disc angular mo- 
mentum vector about the orbital angular momentum vector 
of the secondary. In this case consideration of the linearised 
fluid equations in their three-dimensional form gives the re- 
sponse equation for the zero pattern speed forcing potential 
with m = 1 as (Papaloizou & Lin 1995a) 



dg 



iJ 



dr 1 [!^2(i„i^)2_^2] dr 



(10) 



in which we have prescribed the complex function g, defined 
through: 

—irzfl^g — h tpo- 
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The e-term in (|l^) is prescribed in order to transform the 
location of the singularity due the resonant denominator. 
This amounts to introducing a small phase shift which is 
equivalent to ct for a small shear viscosity, v. The remaining 
symbols correspond to 



pz 



(tpo + 2iu>pr zil sin S)dz, 



(11) 



and 

M = 



2 J 

pz dz 



2n + 3' 



We employ the function g to quantify the degree of warp- 
ing manifest in our numerical calculations. This is possible 
since a constant value of g corresponds to a rigid tilt of the 
disc, for which r ■ V = Q. Then for small amplitude warps 
g — V^/V^ = i(/r, where (^(r) exp(i(;/>) is the vertical dis- 
placement (Papaloizou & Lin 1995a, Larwood & Papaloizou 
1997). 

In the limit of slowly precessing small amplitude warps we 
may neglect the contribution to ( pj| ) from the second term 
in parenthesis and make use of the identity 



-H 2 

pz 



dz = T.n~ 



satisfied then the disc would precess differentially. This ef- 
fect amounts to a situation in which we have the differen- 
tial precession of sonically-connected annuli. For a parabolic 
encounter this would result in a net precession of the disc 
with the main contribution coming from the precession of 
an outer annulus of the disc. 

In addition to the above, Papaloizou & Terquem (1995) 
showed that the j/-component of the tidal torque Ty resulting 
from secular terms Vf' oc 2 gives rise to inclination evolution 
according to 

dS _ Ty 

di ^ 1^' 

Ty is identically zero unless we can introduce a phase shift 
due to viscosity. Since ujp ~ T^/Jo, we notice that the rela- 
tive importance of inclination evolution compared with pre- 
cessional effects in distant encounters is measured by the 
ratio 



Currently, values of a quoted for protostellar discs are typ- 
ically less than ~ 10~^. Therefore we should expect that as 
long as the disc response is linear then precessional effects 
are more significant in determining the final orientation of 
the disc than inclination changes. 



derived by integrating the equation of vertical hydrostatic 
equilibrium (^. Then to determine g we need to integrate 
( p^ over the disc twice. Carrying out this operation assum- 
ing K — Q and a , we find the total range in the vertical 
displacement A{(/r) to be 



2n + 3 f Ro 
5 — n \ q 



— I sin 2(5. 



(12) 



If shear viscosity is a small effect then bending waves prop- 
agate acoustically and should damp on the sound-crossing 
timescale of the disc. This can provide a test of the im- 
portance of any unquantified numerical viscous effects that 
might be introduced due to the vertical shearing motions 
present in warps (see below). 

We note also that in the above calculation there is an im- 
plicit assumption that the sound-crossing time in the disc is 
small enough that the warp is near its steady-state value af- 
ter the encounter. This is equivalent to the requirement that 
the interval over which the interaction occurs is not less than 
the sound-crossing time for the disc. If this condition is sat- 
isfied then it follows that the disc will precess approximately 
as a rigid body (Papaloizou & Terquem 1995) with a small 
frequency cup: 



LUp sin S ■ 



J So ^^r^ sin (fidrdip 
/ Eor3f7drd<^ ' 



where the integrals are to be taken over the area of the disc 
midplane. Evaluating from r = to r = i?o we find 

3 

cos 5. (13) 



Ljp 3 7 — 2n 



4 5 



n \ q 



At large distances of pericentre global rigid body preces- 
sion is possible since the interaction interval is ~ 2n/Qo- 
Larwood et al. (1996) found for circular orbit binaries that 
when the condition for rigid body precession was marginally 



6 NUMERICAL SIMULATIONS 

The numerical method we employ to study the dynamics of 
viscous accretion discs is a version of SPH (see Monaghan 
1992, and references therein) developed by Nelson & Pa- 
paloizou (1993, 1994). The formulation uses a spatially vari- 
able smoothing length, associated with each particle, defined 
in such a way as to ensure accurate energy conservation. 
Each particle's smoothing length is defined to be half of the 
mean distance of the six most distant nearest neighbouring 
particles chosen from a list of forty-five members at each 
time-step. 

In order to stabilize the calculations in the presence of shocks 
an artificial viscous pressure term is included, according to 
the prescription of Monaghan & Gingold (1983). Although 
intended to prevent particle penetration while giving posi- 
tive definite dissipation, the practical implimentation of the 
artificial viscosity introduces a shear viscosity which results 
in disc spreading much as in the standard theory of accretion 
discs (Lynden-Bell & Pringle 1974) . 

These standard accretion disc models were developed to 
study tidal interaction and have been used in previously 
published work (Larwood et al. 1996, Larwood & Papaloizou 
1997). Similar models have been used by Artymowicz & 
Lubow (1994). 



6.1 Basic numerical model 

The SPH disc models were constructed by generating 
random position data for 17500 particles, giving a zero- 
thickness disc of unit radius centred on the origin of the com- 
putational coordinate system (being coincident with the pri- 
mary). Vertical positions were generated in a similar fashion 
by equally seperating this initial set of particles into seven 
uniformly-spaced planes of approximately constant surface 
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density, having total vertical extent Ho- The particles were 
then allocated velocities according to a softened Keplerian 
potential: ^ = (r^ + fc^)"^/*, the softening length b = 0.2 
chosen to avoid the singularity of a Keplerian potential that 
occurs near the origin. It is then only required to determine 
the polytropic constant, K, from the equation of vertical 
hydrostatic equilibrium. We choose K to give a total verti- 
cal semi-thickness Ho at the initial outer radius, given that 
we employ a polytrope of index n — 1.5 in our numerical 
models (which ensures u « constant, Larwood et al. 1996). 
Units are such that G = Mp — 1, giving a time unit of Sl^^. 
The affect of numerical viscosity is to produce an effective 
shear viscosity which we calibrate against a standard con- 
stant viscosity model. The results of our calibration exercise 
imply 



-2/3 



Before introducing the secondary the disc models were al- 
lowed to relax for a few complete rotations solely in the 
gravitational field of the primary and under SPH pressure 
and viscous forces. This is required to establish an approxi- 
mate equilibrium state of the system. Thick disc models with 
H/r = 0.15 (relevant to protostellar discs) as well as thin 
disc models with H/r — 0.05 (chosen to test discs with a sig- 
nificantly smaller sound speed) are considered. Due to the 
different thicknesses and the nature of the SPH numerical 
viscosity the thick and thin disc models also represent dif- 
ferent viscosity regimes (Larwood & Papaloizou 1997), with 
u ~ 10"'^ in the thick disc case and u ~ 10^^ for the thin 
disc models. As a result the outer radius of the disc relaxes 
to different values in each case. For the thin disc models we 
take -Ro = 1.2 and for the thick disc models we use Ro = 1.4. 
The secondary is introduced at a distance of 10 units from 
the primary, corresponding to a time ~ —15, in prograde ro- 
tation with the disc. The secondary's position is given from 
the computed time at each SPH integration step by the ex- 
pressions derived in Section 2.2. The potential due to the 
secondary is then calculated from equation (^). These ex- 
pressions for the secondary are used in their softened form. 
Parameters considered are for secondaries with unit mass 
ratio (Ma = 1) and distances of closest approach q — 2.4, 
3.6, 4.8 and 6.0. Both coplanar and non-coplanar cases with 
(5 = 15°, 30° and 45° are tested. 



6.3 Determining changes to the disc orientation 

As in previous work, we introduce angles t and H to quantify 
respectively the effects of warping and of precession in our 
numerical data. 

J A ■ Jn 



UaIIJi 



and 



cosH 



(Jo X Jd) ■ u 
|JoxJd||u| 



Jd is the disc angular momentum calculated as the sum 
over all disc particle angular momenta, J a is the angular 
momentum within an annulus. The arbitrary reference vec- 
tor u lies in the orbital plane of the secondary. For conve- 
nience we choose u such that H takes the initial value of 7t/2 
radians, its subsequent evolution directly measuring the disc 
precession. 

Retrograde precession of Jd about Jo implies that the an- 
gle n decreases linearly with time, descirbing an arc in the 
orbital plane of the secondary. Applying the pericentre pre- 
cession frequency for one outer-disc orbital period, we esti- 
mate the small change in the disc precession angle expected 
to occur in a 71 = 1.5 polytrope as AH: 

127t f RoV ^ 
An = — 1 cos 0. 



7 

During the parabolic fly-by of a secondary body, the linear 
response of a fluid disc is to precess approximately as a rigid 
body (see below) . The small amount of precession results in 
a shift of the disc midplane away from its initial configu- 
ration, providing an apparent relative tilt with respect to 
the initial midplane, whilst maintaining the angle between 
the disc and orbital angular momentum vectors. For a small 
change in the precession angle AH, producing a small tilt r, 
we expect 

r = |An| sin (5. 

Typical model parameters are fi = 1, q — 4Ro and S — 30°, 
giving a tilt ~ 2°, being much larger than the expected 
change to the inclination S. 

We note that for small values, the total range in t (being 
the angle between the angular momentum vector of the disc 
material contained within a specified cylindrical annulus and 
the total angular momentum vector of the disc) is equivalent 
to A(C/r). 



6.2 Calculation of the disc surface density 

We obtain an azimuthally averaged disc surface density pro- 
file by considering the particle distribution over 140 concen- 
tric circular annuli of width 0.1, taken from radius r = 0.1 
to 1.5 in the midplane of the disc. The surface density at 
the location of the centre of each radial bin is ascribed the 
mean value determined by projecting particles onto the mid- 
plane, summing the number in each annulus and dividing 
this number by the area of the annulus. Surface densities 
thus obtained are expressed in units of the disc's mean sur- 
face density at initialisation. 



7 NUMERICAL RESULTS 

All our disc models lost angular momentum to the sec- 
ondary, resulting in a negative sign in the net angular mo- 
mentum exchanges. Our thick disc models, having the larger 
outer radius, showed the largest loss of angular momentum; 
our thin disc models, having a smaller outer radius, lost 
less angular momentum. Non-axisymmetric waves with az- 
imuthal mode number m = 2 were visibly excited in all but 
the most distant encounters, as we would expect given that 
the m — 2 inner Lindblad resonance with pattern speed 2uJo 
lies far from the disc edge when Ro/q is sufficiently small. 
Additionally, in non-coplanar configurations, we found that 
the disc became warped and showed precessional effects and 
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Table 1. The change in the angu- 
lar momentum content of the disc for 
model parameters used in coplanar 
encounters. 



model 


H/r 


9 


AJ 


1 


0.15 


2.4 


-1.6x10-1 


2 


0.15 


3.6 


-3.0x10-2 


3 


0.15 


4.8 


-2.5x10-3 


4 


0.15 


6.0 


-7.6x10-* 


5 


0.05 


2.4 


-1.4x10-1 


6 


0.05 


3.6 


-1.1x10-2 


7 


0.05 


4.8 


-3.8x10-* 


8 


0.05 


6.0 


-1.9x10-^ 



inclination changes. The affect on the response of introduc- 
ing a finite orbital inclination was generally to diminish the 
effectiveness of the secondary's tide, although there were 
some indications that high inclination encounters could en- 
hance angular momentum transfer. 

7.1 Coplanar encounters 

7.1.1 Angular momentum exchange 

The parameters for these models are given in Table 
(0), along with the corresponding net angular momentum 
changes found in the disc models. Figure (0) shows a par- 
ticle projection plot for the disc particles in model 3, taken 
at a time for which the secondary is close to pericentre. 
The angular separation between the longitude of the sec- 
ondary and the major axis of the disc's elliptical envelope is 
~ 7t/4 radians. If we identify this with the dynamical phase 
shift, rj, then we would have a much larger phase shift than 
would be expected from considering dissipative effects alone. 
Hence the disc viscosity would not be expected to make a 
significant contribution to the process of angular momentum 
exchange in these models. In fact our results indicate that 
the difference in the magnitudes of the angular momentum 
exchanges between the thick and the thin disc models can 
be explained solely in terms of the different radial extents 
associated with the two types of model. 
The magnitude of r] at pericentre is determined by the angu- 
lar velocity of the secondary in the frame rotating with the 
outer disc, the former varying in time. As a result models 
which have smaller (larger) g-values show smaller (larger) 
?7- values. For the range of g-values considered here we shall 
take the approximate median value = 7t/4 in evaluating 
our analytically derived expressions, noting that this gives 
maximal angular momentum transfer. 

We present plots of the total angular momentum of the disc 
versus time for models 5 to 8 in Figure (^). Models 5 and 6 
show an impulsive form for the angular momentum exchange 
which takes place over ~ 27tfl^i time units about pericen- 
tre passage. This simple picture shows signs of breakdown 
in models 7 and 8, which indicate in addition a significant 
return of angular momentum from the secondary to the disc 
occuring after the initial interaction described above. This 
reversal in the direction of the angular momentum transfer 
can be understood in terms of the passage of the outer disc 
through the point where the m = 2 elliptical pattern and 



Figure 1. Particle positions projected onto the {x,y) plane for 
model 3, given at a time for which the secondary is close to peri- 
centre. Each point represents a particle position, the secondary's 
position is denoted by an open circle and its projected path by a 
dotted line. 



the longitude of the secondary are 90° out of phase (giving 
zero net torque at that instant). We can also say that we 
expect the return of angular momentum to be more effec- 
tive for parabolae with larger g-values, since in those cases 
Ro/D varies more slowly in time. Indeed, as the distance 
of closest approach becomes sufficiently large, the net an- 
gular momentum exchange should tend to zero as the disc 
responds adiabatically to the perturbation. In the case of a 
viscous disc the angular momentum exchange would tend to 
the small but finite value expected from viscous dissipation 
alone. 

Our analytical expression for A J does not include the af- 
fect of the return of angular momentum to the disc. We 
compare values derived from equation (|^) with the angular 
momentum exchanges we determine in our numerical mod- 
els, considering only the part of the interaction in which the 
disc initially loses angular momentum (see Figure |^. We 
note that an additional return of angular momentum to the 
disc will modify the magnitude of the resultant exchange by 
a factor which depends on q. 

As we should expect, our analytically derived expression 
gives the best agreement with our numerical results for the 
most distant encounters. However the expression slightly 
over-estimates these values, which is probably a consequence 
of choosing a phase shift that gives maximal angular mo- 
mentum transfer. The numerical values for AJ are under- 
estimated by the predicted values as the encounters become 
closer. This probably owes less to the exclusion of m = 1 
and m — 3 modes from our calculations than to the onset 
of non-linearity in the disc response, which in some cases 
has been found to double the angular momentum transfer 
seen in similar but purely linear calculations (Korykansky 
& Papaloizou 1995). 
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Figure 2. The angular momentum of the disc, in units of its initial value, versus time for models 5-8. 



Figure 3. The magnitude of the angular momentum initially lost 
by the disc for models 1-8. The thick disc data is denoted by open 
squares, the thin disc data by open triangles and the analytical 
result by a solid line. 

7.1.2 The redistribution of disc matter 

In all cases with q ^ 4_Ro, the surface density profile at 
equal times (examined at a total time elapsed of about 30 
units) was found to be barely distinguishable from that of 
an unperturbed model. For q ~ 3Ro there were indications 
that tidal truncation of the disc had occured in much the 



Figure 4. Surface density profiles for model 2, denoted by solid 
circles, and an unperturbed model, denoted by open circles. Each 
is taken for a total run time of about 30 units. 

same way as would be expected for a system involving a 
bound secondary (Lin & Papaloizou 1979a). This case car- 
ries a weakly non-linear component to the response with the 
negative angular momentum fiux into the disc being suffi- 
ciently large so as to result in local non-linear dissipation 
near the disc edge. This behaviour of the system confirms 
that the break-down of linearity occurs for q < 4i?o (cf. Hall 
et al. 1996). The affect of tidal truncation on the disc sur- 
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face density profile is represented in Figure (j^. The middle 
part of the perturbed disc was sculpted into a more com- 
pact and uniformly distributed structure and the outer disc 
was abruptly cut-off. This feature will only be subject to 
significant subsequent evolution on the viscous timescale. 
In models 1 and 5 with q ~ 2Ro the disc response was 
strongly non-linear, the surface density structure being 
severely disrupted and truncated. As indicated in Figure (H) 
a material arm of several disc radii in extent is ejected from 
the outer-edge of the disc and the secondary appears to cap- 
ture some of the disc particles. Much of the remaining disc 
material is re-formed into a dense ring. The formation of a 
ring in both of these models is accompanied by a suppression 
of central mass accretion. This indicates that the ring forms 
by the reflection of an inwardly propagating trailing wave, 
carrying a negative angular momentum flux, producing an 
outwardly propagating leading wave carrying negative an- 
gular momentum from the inner region of the disc. The re- 
flected wave then encounters the rear portion of its inwardly 
propagating counterpart and dissipation occurs. The net ef- 
fect is to clear an inner annulus of material, form a ring and 
suppress central mass accretion. 

The production of circumprimary rings in close binary en- 
counters has been observed before in numerical simulations 
(Clarke & Pringle 1993). It is understandable that the for- 
mation of a ring is only found in very strong encounters 
since it is in these cases that the m = 2 spiral density waves, 
launched at the outer disc, are strong enough to propagate 
to sufficiently small radii that their radial wavelength falls 
below the resolution limit of any code that is used to model 
the disc matter. We find that the resolution at the estimated 
radius of refiection is marginal given the predicted radial 
wavelength (given by equation m; hence we conclude that 
ring formation is observed in these simulations but that it 
is possibly explicable as a numerical effect. 
Accompanying these phenomena we also find that the enve- 
lope of the disc takes on an elliptical shape, however unlike 
the m = 2 bar that is generated prior to strong interaction, 
this form is not centred on the primary. This feature does 
not disperse on a dynamical timescale which implies that 
the disc has become genuinely elliptical. The implication is 
that the disc has taken on an asymmetric structure that will 
only disperse on the viscous timescale. In Figure ^ we plot 
the side view particle projections corresponding to Figure 

(!)■ 

In models 1 and 5 the secondary removed ~ 5% of the initial 
disc mass (i.e. about 10'^ particles), most of which was cap- 
tured. The captured material seemed to consistently mani- 
fest as a dense clump of eccentrically orbiting material (see 
Figure |^. Further investigation of this phenomenon shall 
be the subject of future studies. In model 2 with q ~ 2.6Ro 
only ~ 10 particles were lost. In models 3, 4, 7 and 8, for 
which q ^ 3-Ro, no material was lost from the disc. 



Figure 7. The angular momentum of the disc, in units of its 
initial value, versus time for models 3, 3a, 3b and 3c; represented 
by squares, crosses, diamonds and pluses respectively. 



Figure 8. The angular momentum of the disc, in units of its 
initial value, versus time for models 8 and 8c; represented by 
squares and crosses respectively. 



7.2 Non— coplanar models 

7.2.1 Angular momentum exchange 

The data for our non-coplanar models is given in Table (j^. 
In Figure (^) we show the angular momentum evolution of 
the disc for models 3, 3a, 36 and 3c. As the orbital inclina- 
tion angle is increased the magnitude of the initial transfer 



Figure 9. The magnitude of the angular momentum initially 
lost by the disc for non-coplanar models. The thick disc data 
is denoted by squares, the thin disc data by triangles and the 
analytical result by a solid line. 



© 0000 RAS, MNRAS 000, 



The tidal disruption of protoplanetary discs 



Figure 5. Positional (x, y) data for model 1 taken at a time t = +13. The same data is displayed at two different scales. 



Figure 6. Positional {x, z) and (y, z) data for model 1 taken at a time t = +13. Only particles with r < 1.2 are plotted. 
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Figure 10. Surface density profiles for model 2, denoted by solid 
circles, and model 2c, denoted by open circles. Each is taken for 
a total run time of about 30 units. 

of angular momentum from the disc to the secondary de- 
creases. Hence inclining the orbit of the secondary results in 
a weaker interaction. In addition the magnitude of the an- 
gular momentum returned from the secondary to the disc, 
compared with the net exchange in each case, is also re- 
duced. This is because the mean distance of the secondary 
from the disc edge is larger at higher inclinations. We illus- 
trate an extreme in this behaviour by examining the angular 
momentum evolution for models 8 and 8c in Figure . Al- 
though the initial transfer is largest in the coplanar model 
the net transfer is larger for the non-coplanar model. 
As for the coplanar case we compare the inferred A J values, 
determined for the initial phase of the angular momentum 
exchange, with our expression derived from linear theory and 
modified for inclined orbits of the type considered here. We 
present these results in Figure (^). The general behaviour 
of the numerical values compared with the analytical result 
is similar to that found in the coplanar case. In addition it 
appears that low inclination models give better agreement 
with the analysis than the high inclination models. This is 
probably due to the enhancement of the angular momen- 
tum flux through the non-coplanar m — 1 contribution to 
the response, which is most important at high inclinations 
(Papaloizou & Terquem 1995). We note also that the thin 
disc models appear to show larger values than the thick disc 
cases, this is probably due to a greater tendency to non- 
linearity in the response which could be more significant in 
the presence of the vertical shearing motions found in warps. 

7.2.2 The redistribution of disc matter 

Tidal truncation of the model discs was found to operate for 
non-coplanar models as for the coplanar cases, but occurring 
at a reduced level for larger orbital inclinations. In Figure 
( p^ we compare the surface density profiles for models 2 
and 2c, examined at a total time elapsed of 30 units from 
the time of the introduction of the secondary (cf. Figure 
The surface density profile for model 2c appears to be 
intermediate between the coplanar case and the unperturbed 
model with the outer disc not as strongly cut-off as in the 
coplanar case but more so than in the unperturbed case 
(which is not cut-off at all) . 



Figure 11. Precession angle data versus time for model 3c. 

Mass transfer streams and spiral arms were generated as in 
the coplanar models except that they could become warped 
and even wrapped (cf. Clarke & Pringle 1993) about the disc, 
for the closest encounters, as shown in Figure (|l2|). Mass loss 
from the disc was also of a similar order to that observed in 
the coplanar models. 

7.2.3 Precession, warping and inclination changes 

All our non-coplanar models showed a small amount of pre- 
cession. In Figure (|l^) we give the precession angle evolution 
for model 3c. As for the torque which gives rise to angular 
momentum exchange the torque giving precession of the disc 
is found to assert itself essentially as an impulse at pericen- 
tre. The small amount of precession manifest in all these 
models was found to give a small relative tilt r ~ 1° — 6°, 
even at the largest g- value, in reasonable agreement with the 
expression derived above. We give particle projection plots 
of model 6c in Figure (p^. It is clear that the midplane of 
the disc has become shifted with respect to its initial orien- 
tation. 

As well as showing precessional behaviour the model discs 
also became warped. We compare the size of the warps found 
in each model at a time for which the companion has re- 
ceeded to r = 10, with the equation (|l^), in Figure (|l^. 
The expression derived from linear theory seems reliable for 
distant encounters in the thick disc models. In the thick 
disc cases with q < 4i?o the magnitude of the warp is over- 
estimated, consistent with non-linear damping of strong ver- 
tical shearing motions. 

The size of the warps present in the thin disc models seem 
to obey a proportionality relation with the tidal strength of 
the encounter but are over-estimated by a factor of ~ 10. 
The three times larger sound-crossing time (= J c~^dr) in 
the thin disc models prevents saturation of the warps on 
the encounter timescale. In Figure (^^ we give the warp 
evolution for model 3c. A decay timescale of approximately 
60 units is inferred, being about 5 sound-crossing times. A 
similar result was found to hold for the thin disc models. 
This is consistent with bending modes that propagate and 
damp as acoustic waves (Papaloizou & Lin 1995a). 
Inclination changes to the disc were typically small in cases 
such that q > 3Ro, for q ~ 3i?o and q ~ 2Ro the inclination 
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Figure 12. Positional data for model Ic at a time t = +13. 



Figure 13. Positional data for model 6c, the top frame gives the unperturbed particle configuration and the bottom frame gives the 
positions at a time t = +28 after pericentre passage. 
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Table 2. Data for model parameters used in non-coplanar encounters. &6 denotes 
the change to the orbital inclination of the disc, all other symbols are defined 
in the text. Angles are expressed in radians except for 5 and the quantities in 
parentheses which are expressed in degrees. 



model 


ri /V 


q 


c 



A 7 








ATT 
Ziil 


T 






la 


0.15 


2.4 


15 


-1.5x10" 


1 


1.7x10" 


2 


-0.19 


4.9x10" 


2 


(2.8) 


lb 


0.15 


2.4 


30 


-1.3x10- 


1 


2.8x10" 


2 


-0.17 


8.5x10" 


2 


(4.9) 


Ic 


0.15 


2.4 


45 


-9.7x 10" 


2 


2.7x 10" 


2 


-0.15 


1.1 x 10" 


1 


(6.1) 


2a 


0.15 


3.6 


15 


-2.0xl0~^ 


2 


4.4x 10" 


3 


-0.11 


2.8x 10" 


2 


(1.6) 


ZD 


0.15 


3.6 


30 


-1.9 X 10~ 


2 


6.1 x 10" 


3 


-0.10 


5. Ox 10" 


2 


(2.9) 


2c 


0.15 


3.6 


45 


-8.7x10" 


3 


5.0x10" 


3 


-0.08 


5.7x10" 


2 


(3.2) 


3a 


0.15 


4.8 


15 


-2.1x10" 


3 


3.6x10" 


4 


-0.08 


2.1x10" 


2 


(1.2) 


3b 


0.15 


4.8 


30 


-1.6x10" 


3 


4.9x10" 


4 


-0.07 


3.5x10" 


2 


(2.0) 


3c 


0.15 


4.8 


45 


-1.1x10" 


3 


2.7x10" 


4 


-0.05 


3.5x10" 


2 


(2.0) 


4a 


0.15 


6.0 


15 


-2.0x10" 


4 


0.6x10" 


4 


-0.06 


1.6x10" 


2 


(0.9) 


4b 


0.15 


6.0 


30 


-1.9x10" 


4 


0.7x10" 


4 


-0.05 


2.5x10" 


2 


(1.4) 


4c 


0.15 


6.0 


45 


-1.5x10" 


4 


0.3x10" 


4 


-0.04 


2.8x10" 


2 


(1.6) 


5c 


0.05 


2.4 


45 


-7.2x10" 


2 


2.7x10" 


2 


-0.15 


1.1x10" 


1 


(6.1) 


6c 


0.05 


3.6 


45 


-4.1x10" 


3 


1.4x10" 


3 


-0.07 


4.9x10" 


2 


(2.8) 


7c 


0.05 


4.8 


45 


-2.4x10" 


4 


0.2x10" 


3 


-0.05 


3.5x10" 


2 


(2.0) 


8c 


0.05 


6.0 


45 


-5.3x10" 


5 


0.3x10" 


3 


-0.04 


2.8x10" 


2 


(1.6) 



Figure 14. The magnitude of the range in the vertical elevation 
at a time of about 30 units. The thick disc data for models with 
q ^ 4i?o is denoted by open squares, thick disc data for models 
with q ~ 3Ro by open circles, the thin disc data by asterisks and 
the analytical result for the thick disc models by a solid line. 



Figure 15. The warp evolution shown as the inclination versus 
radius for approximately equal time intervals. Squares represent 
data at time t = -|-31, pluses represent data at time t = +39, 
triangles represent data at time t = +50, crosses represent data 
at time t = +61 and diamonds represent data at time t = +70. 



changes were more significant at ~ 1% and ~ 10% respec- 
tively. The net change to the inclination was found to be 
much smaller than the tilt due to precession in all cases. 

8 DISCUSSION 
8.1 Summary 

In this paper we have presented the results of analytical 
calculations and hydrodynamic simulations of the tidal in- 
teraction of a point mass encountering a viscous accretion 
disc on a parabolic prograde fly-by. We have considered the 
interaction for the case when the orbital plane of the sec- 
ondary coincides with the midplane of the disc as well as 
the case when the two planes are inclined. 



We have shown that for coplanar systems the most impor- 
tant component of the perturbing potential that gives a net 
tidal effect is the m — 2 bar mode applied at pericentre. 
Further to this we have demonstrated that the affect of in- 
troducing an anomolous shear viscosity of the magnitude 
thought to be relevant to protostellar accretion discs is only 
significant for very distant encounters. The exchange of an- 
gular momentum occuring between the secondary and the 
accretion disc is calculated from simplified inviscid consid- 
erations by modelling the interaction as a non-resonant im- 
pulse occurring at pericentre. A net torque is found to act 
on the disc due to the generation of a phase lag in the den- 
sity response. To calculate the disc response we apply the 
pericentre torque over a finite interaction interval for a fixed 
value of the phase shift. 

In our non-coplanar models we have shown in addition that 
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a m = 1 bending mode is excited whose action is to warp 
the disc, cause it to precess away from its original plane and 
give a small change to the inclination. The precession of the 
disc produces a small tilt of the disc midplane relative to 
its initial configuration, this effect is always more important 
than the inclination change in determining the final orien- 
tation of the disc. The size of the tilt is significant even for 
the weakest encounters considered, namely those with the 
distance of closest approach of order five times the outer 
radius of the disc, for which the exchange of angular mo- 
mentum is a very small fraction of the total disc content. 
There is also evidence to suggest that the m — 1 component 
of the perturbing potential can become as important as the 
m — 2 component in the angular momentum exchange at 
high inclinations. 



If we then consider the simplified case with the mass of the 
disc distributed in a thin ring of radius _Ro, being in Keple- 
rian rotation about the primary alone, then we find: 

Then considering further extreme parameters for angular 
momentum exchange, such as /i ~ 1, e 1 and q ~ 2Ro we 
find that 6e ~ AJ. Thus if we take the values for AJ present 
in our models we would expect that eccentricity changes to 
the secondary are generally less significant. The implication 
of this is that an eccentric binary may remove a significant 
fraction of the disc's angular momentum without becoming 
unbound. 



8.2 Consequences for protostellar accretion discs 

8.2.1 Silhouette discs in the Orion Nebula 

Clarke & Pringle (1991) estimated the encounter rate, F, for 
parabolic point masses encountering discs. For the densest 
parts of the Trapezium cluster they estimated F ~ O.lMyr"^ 
for discs with Ro ~ lOOAU. They argued that since disc 
lifetimes are ~ IMyr then disc penetrating encounters are 
rare. We note that approximately F cx q^. The implication 
of this is that encounters with q ~ 3i?o are the rule rather 
than the exception. 

In our numerical models encounters with q ~ 3Ro produce 
disc tilts ~ 3° and show tidal truncation such that the discs 
become more compact with strongly cut-off edges. Directly 
imaged silhouette discs in the Trapezium cluster are found 
to exhibit a similar property in their surface density profiles 
(McCaughrean & O'Dell 1996). Also the largest silhouette 
disc {Orion 114-436; with Ro ~ 500AU) which is almost 
edge-on to the line of sight, shows a clear radial asymme- 
try and vertical distortion. But note that other environmen- 
tal influences cannot be ruled out as explainations for these 
features of the observations (O'Dell, Wen & Hu 1993, Mc- 
Caughrean & O'DeU 1996). 



8.2.2 Discs in eccentric binary systems 

Korykansky & Papaloizou (1995) noted that it is possible 
to model the accretion disc response in an eccentric binary 
system with a succession of parabolic pericentre passages. 
In this picture the accretion disc suffers a series of impacts 
in which it loses angular momentum to the secondary, thus 
reducing its lifetime. The secondary would increase its ec- 
centricity with every impact. 

In coordinates based on the primary, a secondary with spe- 
cific angular momentum L and orbital eccentricity e has a 
pericentre distance q: 

^ " 1 + e G(Mp + Ms) ■ 

Therefore at pericentre we can write the change in the ec- 
centricity 5e in terms of the change to the specific angular 
momentum content of the disc 6 J: 

6e = -2(l + e)^. 



8.3 Consequences for planetary systems 

8.3.1 The Solar system 

If planets form near the midplane of an accretion disc and if 
protoplanetary discs can become tilted by some mechanism 
without exerting a significant torque upon the central star, 
then it is reasonable to suppose that planetary systems are 
able to form with non-zero inclinations to the stellar equator. 
We note that it is trivial to show from standard astronomical 
data that the invariable plane of the Solar system is tilted 
by ~ 6° with respect to the Solar equator. 
The largest tilt angles we found in our simulations were ~ 6° 
for q ~ 2Ro. It is therefore plausible that a unique close 
encounter in the early history of the Solar system could have 
caused the Solar obliquity. If on the other hand we were to 
suppose a number of more distant but coherent encounters 
then we must appeal to the existence of an eccentric Solar 
binary companion which has long since become unbound 
through its tidal interaction with the primitive Solar disc. 
However, our result of the previous Subsection indicates that 
unless it had an eccentricity very close to unity, a Solar 
binary should still be in evidence. 

8.3.2 The Beta Pictoris system 

A gaseous disc in which there already exists a planet lo- 
cated within a tidally cleared gap (Lin & Papaloizou 1979b) 
interacts with that planet through gravitational torques. 
If the outer disc becomes tilted due to a close binary en- 
counter then the planetary orbit would not necessarily tilt 
rigidly with the rest of the disc. A relative inclination be- 
tween the disc and the planetary orbit would be generated. 
We note that the dusty circumstellar disc associated with 
Beta Pictoris shows an inner warp ~ 3° (Burrows, Krist & 
Stapelfeldt 1995). This has been modelled with an inclined 
planet orbiting interior to the disc (Mouillet et al. 1997). 
Asymmetries on the scale of the disc have also been noted 
and a stellar encounter proposed as a possible mechanism 
for their generation (Kalas & Jewitt 1995). 
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